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Abstract 

In this note we propose an exact simulation algorithm for the solution of 



Oh! dXt = dWt + b{Xt)dt, Xo=x, (1) 

, where b is a smooth real function except at point where b(0+) 7^ b{0—). The 

main idea is to sample an exact skeleton of X using an algorithm deduced from 
the convergence of the solutions of the skew perturbed equation 

(N 

> . dXf = dWt + b{Xt^)dt + l3dL'l{X^), Xo^x (2) 

as _ 

■ towards X solution of ([T]| as /? 7^ tends to 0. 

CO ■ In this note, we show that this convergence induces the convergence of exact 

^ - ^ simulation algorithms proposed by the authors in [7] for the solutions of ([2]) 

fT^ ' towards a limit algorithm. Thanks to stability properties of the rejection 

procedures involved as /? tends to 0, we prove that this limit algorithm is an 
exact simulation algorithm for the solution of the limit equation IT}. Numerical 
examples are shown to illustrate the performance of this exact simulation 
algorithm. 
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1. Introduction 

1.1. Motivations and exposition of the problem 

Exact simulation methods for trajectories of one-dimensional SDEs has been a subject of much interest 
in the last years : see for example [S], [3], [3], [H], [T3]. Unlike the classical simulation methods, which all 
involve some kind of discretization error (see for example [Ij for the Euler Scheme), the exact simulation 
methods are constructed in such a way that they do not present any discretization error (under the strong 
hypothesis that the diffusion coefhcient is constant and equal to one). In the last years, the original method 
presented in the fundamental article [5] has been extended to overcome various limitations of the initial 
algorithm ; it has been generalized to include the cases of unbounded drifts ([3], [Sj) and extended to various 
'non classical' type of SDE (j7j). 

In this paper, we present an attempt for the adaptation of the exact simulation methods of [5] to one- 
dimensional SDE that possess a discontinuous drift at point 0. Namely, our object of study is {Xt)t>o 
solution of 

dXt = Wt+ h{Xt)dt, Xo = X, (3) 

where 6 is a smooth real function except possibly at point where ^(0-1-) ^ 6(0—). 

The simplest case of a process solution of an equation of type ^ is surely the so-called 'Brownian motion 
with two valued drift' solution of 

dXt^Wt + i0olx,>o + Oilx,<o)dt, Xo = x, (4) 

where {0ci,di) £ M^. For a general reference concerning these type of motions, we refer to [TD] p. 440-441 
or [Hj. These motions appear in stochastic control problems (see for example [2], [^) and also theoretical 
studies concerning representations of reflected Brownian motion with drift (see [S] in the case 6*0 = — ^i). 
Even though there exists explicit representation formulae for the densities of such Brownian motions with 
two valued drift in terms of combination of convolution integrals (see [lOj p. 440-441), up to our knowledge 
there is no exact numerical simulation algorithm for such motions available in the literature. The algorithm 
presented in this paper gives an answer to this question. 

1.2. Main ideas of the paper 

In ^7j, the authors manage to adapt the exact simulation methods of [3J to the case of one-dimensional 
SDEs that possess an additional term involving the local time of the unknown process at point 0. Namely, 
the exact simulation methods of [3J are modified in [7i to include the case where (Xf )(>o is the solution of 

dX^ ^Wt + b{X^)dt + l3dL°iX^), Xo = x. (5) 

In this situation 7^ |/?| < 1 (with (3 ^ 0). L^iX) denotes the symmetric local time of X^ in zero at time t, 
and h is still allowed to be discontinuous at 0. 
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The main idea in [7J was to propose an exact rejection simulation algorithm for the solutions of ([5]) 
using as sampling reference measure the law of some drifted skew Brownian motion with prescribed terminal 
distribution and with drift of magnitude 1//3 avoiding the case where /3 = 0, for which we propose a proper 
treatment here. Our contribution in [7] deals mainly on the simulation of such bridges of drifted skew 
Brownian motions using a classical rejection procedure and looking for tractable rejection functions. 

Unfortunately, a direct exact simulation method along the same lines as fj] cannot be properly defined in 
the case where /3 — 0. However, we know from Le Gall in [TT] that solution of (O tends strongly to X 
the solution of (jS]) as /? tends to on each time interval [0,r]. This leads us to examine what happens at 
the level of the algorithms proposed in [7J as /3 tends to 0. 

In fact, we check here by computations that there is indeed a convergence phenomenon at the level 
of rejection functions and rejection sets involved in the exact simulation algorithms given in [7J. This 
convergence gives rise naturally to a nice and implementable limiting algorithm. 

The main problem becomes then to prove rigorously that this limiting algorithm is indeed an exact 
simulation algorithm for the solution of ([3]). In particular, as far as we see, the direct interpretation of 
this limiting algorithm is not clear ; for the time being, we have to confess that we really understand the 
construction of the limiting algorithm exposed in this paper only via the convergence procedure explained 
above. Let us also emphasize that this new algorithm is still a rejection algorithm, and one may naturally ask 
for a direct interpretation of its corresponding reference measure. In Remark 13. II we give an interpretation 
of the reference measure (corresponding to the limit rejection algorithm) in terms of a standard Brownian 
motion conditioned on prescribed laws for its final position and its local time at at time horizon T. 



1.3. Outline of the paper 

The paper is organized as follows. In the preliminary Section 2, we explain the convergence of rejection 
sampling algorithms in a general framework. The result exposed in this section will be used to justify that 
our limiting algorithm is indeed an exact simulation algorithm for the solution of ([3]). The exact simulation 
problem treated here is presented in Section 3, where we explain the manner in which we adapt the exact 
simulation methods of [3J to our situation. Yet, the resulting algorithm adapted from [3J is not directly 
implementable in our context because we have to sample from a complicated reference probability measure 
Z. The sections 4 and 5 are devoted to the interpretation of Z as a limit of some sequence ^Z„^ of better 
known probability measures. Finally in Section 6, we apply the results of the preliminary Section 2 to the 
sequence ^Z„^ . This gives rise to a directly implementable limit algorithm for the exact simulation of a 
skeleton along the reference probability Z. We end up the article with numerical results and illustrative 
examples shown in Section 7. 
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2. Preliminary : convergence of abstract rejection sampling algorithms 

Proposition 2.1. i) Assume that we have a sequence (^„) of probability measures on a measurable space 
{S,S), and ^dom o- probability measure on {S,S), satisfying for any n G N 

d^n 

with £„ > and < /„ < 1. 

Assume that /« — > / as n —> oo point-wise on S . 

Then, (^„) converges towards a probability measure ^ satisfying 



with e = lini„^oo £« • 

ii) Moreover, let {Yk, Ik)k>i be a sequence of i.i.d. random elements taking values in S x {0,1} such 
that Yi ^ S,dom o,nd P[/i = l|Yi = y\ = f{y) for all y (z S. Define r = min(fc > 1 = /fc = 1). Then, 
V{Yredy)^ady). 

Proof. For any A £ S we have £,n{A) = Jj^ fn{z)^dom{dz). By dominated convergence we have 

fn{z)Uom{dz) > / f{z)£,dom{dz). 

A "-^"^ Ja 

Taking A = S, and as £,n{S) = 1 for any n e N, we have 

1 1 



e. 



Js fn{z)Uom{dz) n^oo Jg f {z)^dorn{dz) 

Setting now for any A £ S, ^{A) — ^ f{z)^dom{dz), it is clear that 

VAeS, UA) ^e(^). 

n— f oo 

Then ^ is a probability measure on {S,S). It satisfies © by construction. This proves point i). For the 
proof of point ii), see Proposition 1 in [6]. 

3. Exact sampling algorithm for an SDE with discontinuous drift (inspired by [6]) 
3.1. Assumptions 

The function 6 : M M is bounded with bounded first derivative on ]R*^+ and R*'~ with a possible 
discontinuity at point {0}. We suppose that both limits lim^^o+ ^(-2) = b{0+) and lim^^o- K^) = ^(0^) 
exist and are finite. The value 6(0) of the function 6 at is of no importance and can be fixed arbitrarily to 
some constant (possibly different from either 6(0+) or 6(0—)). 

We introduce the notation 

^ _ 6(0+) - 6(0-) 
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3.2. Change of probability 

Let < T < oo. Denote C — C([0, T] , R) the set of continuous mappings from [0, T] to R and C the Borel 
cr-field on C induced by the supreme norm. 

Let P be a probabihty measure on (C, C) and W a Brownian motion under P together with its completed 
natural filtration {J^t)t>o- We will denote P^ = P(- 1 Wq — x). When necessary we will denote by w = 
(a;t)o<t<T the coordinate process. 

We consider the following SDE 

dXt = dWt + h{Xt)dt, Xo = X. (8) 

Our objective is to sample along Xt- 

Let us define on (C, C) the probability measure W by 



^ = exp{-^ -biX,)dW,-'-l l^iX,)dt}. 



Under W the process X is a Brownian motion and we have, 

i-T 



^=exp{£B(W-i£p(X.).t}. 



Thus for any bounded measurable functional F : {C,C) R we have, 



T 1 '■^ 



E^[F{X)]^E^[F{X)e^p{ I b{Xt)dXt~-j^ ¥{Xt)dt)]. (9) 







We set B{x) = b{y)dy. Using the symmetric Ito-Tanaka formula (see Exercise VI-1-25 in [T5]), and 
the Occupation times formula ([13J) we get 



T -, pT 



B{Xt) - B{Xo) = I b{Xt)dXt + \l b'{Xt)dt + ^0+) ^ ) j,o^(x). 







Thus ^ becomes, 

E^[F{X)]^E^[F{X)exp{B{XT)~B{x)-eL''r{X)- f (l>{Xt)dt}], 

Jo 

where we have set 

P{x)+b'{x) 

nx) = . 

Setting now 

(j){x) = 4>{x) — inf 4>{x), 

a: (EM 

we finally get that for any bounded and measurable functional F : (C, C) ^ R we have, 

E^[F(X)]cxE^[i^(X)exp{B(XT)-B(.T)-^?LO(X)}exp{~ / ^{Xt)dt]]. 

Jo 
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Let us now introduce the probability measure Z on (C, C) defined in the foUowing way 

^(c^) (X exp - B{x) - eL^^{X){u)). 

We have 

E^[F(X)] cx Eg [F{X) exp { - f ~^{Xt)dt]] . 

Jo 

In the sequel we note Z the probability measure induced on (C, C) by the law of X under Z. 

Remark 3.1. (Interpretation of the probability Z) 

Recall that under W the process X is a Brownian motion and that, by definition, 

^(^) (X exp {B{Xt{u:)) - B{x) - 

In particular, under the probability Z, X is a Brownian motion conditioned on [XttL^) ~ h{y,£)dyM with 

h{y,i)dyd£ cx exp {B{y) - B{x) - 91) W {Xt G dy, L% e di) . 

This makes it difficult to sample exactly Xt under Z for t e (0,r). 

3.3. Exact simulation algorithm for the solution of ^ starting from x 

Let us denote by K an upper bound for 4'{x). Following the spirit of [3j we can thus sample from Xt 
using the following algorithm. 

EXACT SIMULATION ALGORITHM FOR THE SOLUTION OF (H]) starting from x. 

1. Simulate a Poisson Point Process with unit density on [0,T] x [0,X]. The result is a 

RANDOM NUMBER N OF POINTS OF COORDINATES (ii , Zi ),..., (t at , Zat) . 

2. Simulate a skeleton {ujh , • ■ • , wt„ , ujt) where a; ^ Z. 

3. If Vi e {1, . . . , N} ^(w*.) < accept the skeleton. Else return to step 1. 



This algorithm produces an exact sampling of Xt under P: it is the final instance lot of an accepted 
skeleton. 

The main issue in the above algorithm is to sample a skeleton of the canonical process under Z (Step 2). 
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Remark 3.2. (Other exact simulation algorithms) 

Other probability changes might be performed in order to try to tackle the exact simulation problem 
presented in the introduction. For example (though we will not prove it here) it is possible to swap to a 
probability measure § under which X has the law of a some Brownian motion with a symmetric two valued 
drift (solution of equation Q in the case where = — ^'i) with some prescribed terminal law. Even though 
the density probability distribution of such bridges may be explicitly computed, it seems difficult to find 
tractable general rejection bounds for these laws. 

4. Recalls on the skew Brownian motion with drift 

In this section, we recall some basic facts concerning the skew Brownian motion with constant drift. 
Although these facts seem at first quite far away from our purpose, they will be used in the sequel in order 
to justify that the limit rejection algorithm presented in Section 6 returns an exact sampling under Z. At 
the end of this section, we give an algorithm for the simulation of bridges of SBM with constant drift that 
will be used as a basic building block in the sequel. 

4.1. The transition function of the skew Brownian motion with drift 

Let us recall that the Skew Brownian Motion (SBM) with constant drift component ji G M., denoted by 
B^'f^, solves 

dB^^" = dWt + iidt + ^dL'^tiB'^'"). 

This SDE with local time has a unique strong solution as soon as |/3| < 1 (see pTTj). The process B^'^ 
enjoys the homogeneous Markov property. We shall denote by p^''-'-{t, x, y) its transition function. 
Let us introduce the function v^'f^{t, x, y) defined by 

vP'^^{t,x,y) = (l-exp(-^))l,,>o 

(10) 

+ (1 + Sgn(2/)/3) exp(-^l,,>o) [l - /3^V2^exp{^^^^i^±^}7V-(^ii^:H^)] , 

where N^iy) ^ e-'^/^dz. 

With this notation we can rewrite the expression of p^'^(t, x, y) given in [7]. 

Proposition 4.1. We have for all t > Q, all x, y G K, 

pP-^{t,x,y)=p'^-^{t,x,yy'>^{t,x,y). (11) 



Proof. See f7] (Proposition 4.7). 
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4.2. Bounds for the transition function of the SBM with drift 

In this paragraph, we give bounds on the transition function of the SBM with drift. These bounds will 
be used in the sequel to find tractable rejection bounds for our algorithm. 
Let us set a = max(ii^, and 

z) = 1 - /3,y2^exp(ii±|^)A.^(^). (12) 

We also set 

, f 2a if /3u > 

cf ;r = < (13) 

[ 2a7^'''(i, |a;|) if /3/Lt < 0. 

We have the following result. 
Lemma 4.1. Let {(3,^) e (-1,1) x R. We have 

v^'^{t,x,y)<4f, Vx,yeR. (14) 

Proof. Equation ([14]) comes from p4|) and the fact that, if Pfi > 0, we have p^'^^{t, x, y) < 2ap^''^{t, x, y) 
for aU x,y eR, and if < 0, we have p^''^{t,x,y) < 2aj^''^{t, \x\)p'^-''^{t, x,y) for aU x,y eR (see, in f7], 
Lemma 5.3 and its proof). 

We also have the following lemma. 

Lemma 4.2. Let {l3,fi) e (-1, 1) x M. We have 

vP'^(t,x,y)<c^,;^, Vx,y€R. (15) 

Proof. This comes again from pT|) . together with the fact that p^'^{t,x,y) < Cf''yp'^'^{t,x,y), for all 
a;,?/ G R (see again, in Lemma 5.3, especially the proof of Equation (5.7)). 

Remark 4.1. Note that v^'^{t,x,y) > and 7^'^(i, z) > for any t £ R*'+, x,y,z e R, even for large 
values of (see Remark 4.8 in [3). 

4.3. Sampling bridges of the SBM with drift 

We denote by q'^''^{t,T,a,b,y) the density defined (for i < T) by 

P[Sf^^ e dy I B^''^ = a, B^^^ = 6] = /'"(i, T, a, 6, y)dy. 

The function {t,y) i— > q^''^{t,T,a,b,y) is the transition density function of a bridge of an SBM with drift 
relating points a and b in T unit time. 

As (?°'^(t, T, a, 6, y) = g"^°(t, T, a, 6, y), by Proposition [H] we get, 

0,^f,rr. , ^ 0,0.,^ , y''^{t,a,y)v^^^{T-t,y,b) 

qP'^'{t,T,a,b,y) = q"'''{t,T,a,b,y) j— — . (16) 

yP'^'ll , a, b) 
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Let us set 

if p^>0 

Aa^"fl^^''{t,\a\)jP'^{T -t,\b\) if P^i<0. 



CtT'a,t={ ...... . . . *: (17) 



We have 

gO.O(t,r,a,6,y) t;^*^A.(r, a, &) ^'^■^ 

with 

.'B.^.p, ,_ v^'>^{t,a,y)v^''^{T~t,y,b) 

Ja,b yV) - „/3,;x ' ^'^'^i 

'~^t,T,a,b 

where the superscript *B appears for the word "Bridge". 
Considering (US]), ([H]), (HHD and ([T7D it is clear that 

We thus propose the following rejection algorithm in order to sample along q^''^{t, T, a, b, y)dy. 

Auxiliary Algorithm 1: Sampling along q'^''^{t,T,a,b,y)dy 

1. Sample a Brownian bridge Y along q°'°{t,T,a,b,y). 

2. Evaluate 

3. Draw U ^ U{[0, 1]). If U < f^f''^{Y) accept the proposed value Y. Else return to Step 2. 



Remark 4.2. Note that the quantities v'^''^, 7^^'^, c^f, C^t a &' ^^^^ fTf^^ defined respectively in ([TU)) . p^ . 
(fT5)) (fT7|) . and involved in the above algorithm depend only on ji through the product /3/i. This 
computational fact gives the key ensuring the construction of the limit algorithm by convergence performed 
at the beginning of Section |6l 

5. Convergence of a sequence of probability measures towards Z 

In this section, for any n G N we denote by AT" the solution of 

dAt" = dWt + b{X^)dt + ^L?(A"), Aq" = x. (19) 
For the existence and uniqueness of solutions to (|19p see [11] . 
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5.1. A strong convergence result 

Theorem 5.1. (Le Gall [llj, 1984.) Let X he the solution of ^ and (X") the sequence of solutions of 
(frail. We have for allO<t<T, 

E[ sup \X,-X^\] ^0. 

0<s<t 

Proof. We use the notations of [TT] . Using the Occupation times formula we can rewrite Equation ([19]) as 

dX^ ^dWt+ f Mdy) dmx^), XS = X, 



with Vn{dy) — b{y)dy + j^So{dy). Lemma 2.1 in [TT] asserts that there is for each n € N a function f^^, 
unique up to a multiplicative constant, satisfying fl^^{dy) + (/j/„(y) + fvn{y~))^n{dy) , where the notation 
fl^^{dy) is for the derivative of /^^ in the generalized sense. Lemma 2.1 in [11 also asserts that if we require 
that fjj^ [x] — > 1 as a; — > — oo then, 

ry _ 1-1/n 
f>^Av) = exp ( - 2 / h{z)dz) X ly>o . 

The sequence of functions (/^^ )„ clearly converges point- wise to /(y) — exp ( — 2 b{z)dz^ . By dominated 
convergence we have for all if > that / ^, - f\{y)dy ^ as n oo. Thus Theorem 3.1 in [llj asserts 
that 



with X the solution of 



E[ sup \Xs-X^\] ^0, 

0<s<t 



dXt =dWt+ [ v{dy) dL\{X), Xa = x, 
Jr 



where i^{dy) ~ ^J(y)+f^P^ ~ ~1 /fa)^ '^^ ~ ^{y)dy- That is to say X is the solution of ([8]). 
5.2. The probability measure Z as a limit of a sequence (Z„)„ 
Recall the definition ([7]) of 9. Let us set 

and bn{x) = b{x) — iin- 
From (fini we have, 

XJ}.=x + H/^^'" + M„T + iL!^(X"), 
where the process W^^''^'' defined by 

^■^SD.n ^ ^ br,{X^)dt 
is a Brownian motion under the probability measure W^^'" defined by 



exp{-^ b^{Xl^)dWt~\j^ bl{X^)dt]. 



(21) 
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Let us now set Bn{x) — bn{z)dz. As shown in [7j using a symmetric Ito-Tanaka formula we have for 
any bounded measurable functional F : (C, C) — > M, 

E^[F(X")] - E^sz,,„ [F{X") cxp {B„(X?) - B„{x) - / (/'„(X,")di}] , 

Jo 

where 0„(z) = ^"(^) + ^"(^^ + ^^"^"('^) . Notice that the law of X" under W^-°'" is the one of a SBM with 
drift. 

Let us define the probability measure Z„ on (C, C) by 

^^^H(xexp{i?„(X^H)-i?„(a;)}, (22) 

and Z„ the probability measure induced on (C, C) by the law of X" under Z„. 

The law Z„ can be well described (see Subsection I5.3p . Thus in [7j we managed to sample exactly along 
(|19p using skeletons under Z„ as proposals and exp { — {(j)n{ujt) — inf 0n)c?i} for a rejection rule. 

For our present purpose it turns out that we have the following result. 
Proposition 5.1. We have 

Z„ — ^ Z. 

Proof. Let be F : (C, C) ^ M a bounded measurable functional. We have from (|22p . (piT) . 

E| [F(w)] =E^JF(X")] 
=c„E^sz,,„ [F{X") exp{B„(X?) - B„(a:)}] 

=c„Ej[F(X")exp{i3(X?) - i3(a;)}exp{-^„(X? - x) - bniX^dWt ~ ^ C bl{X^)dt)], 

Jo ^ Jo 

where c„ is the normalizing constant ensuring that Z„ is a probability measure. But 

- /i„(x^ r b^ixj^)dwt -\ f bi{x^)dt 

Jo ^ Jo 

= - pl^Wt - [ finbiXj')dt - -finL'^TiX'') ~ [ biX^)dWt + finWr 
Jo n Jo 

^ ^ ¥{X^)dt+ I ^iJ{X^)dt-lplT 

T 



n Jo 2 Jo 2 



Thus, 



= c„e-^Ei^[F(X")exp{B(X?)-B(x) - ^L?,(X")}exp{- / 5(Xt")dM^t - ^ / P{XJ')dt} 

n Jo 2 Jo 
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Suppose first that z ^ B{z) is bounded. From (|20|) . we see that ^ 6* as n oo. and combining the 
results of Theorem 15 . II and Lebesgue's convergence theorem for stochastic integrals (see for example Theorem 
(2.12) p.l42 in [13]), we have 

hm Eg [F(^)] 

= liiR{cne-^)E^[F{X)exp{B{XT)-B{x)-eL°r{X)}exp{- [ h{Xt)dWt-\f b'{Xt)dt)]. 

When z i-^ B{z) is not bounded, one may use a truncation argument in order to retrieve the above equalitiy. 
By definition of W, 

lim Eg [f^(^)] 

n— !-oo " 

= lim (c„e-^)E^[F(X)exp{B(XT) - B{x) - eL%{X)}] 

n— )-oo 

= hm {cne-^)W^[F{X)] 

n—^oo 

= lim (c„e-4^)E|[F(c^)]. 

For the special case _F = 1, the above computation ensures that \\ni^^^^{cjie 2—) = 1. Thus, for any 
bounded and measurable F : (C, C) ^ R we have 

lim E| [F{uj)]^W^[F{lo)] 

and the proposition is proved. 

5.3. Sampling a skeleton under Z„ 

We have the following proposition. 

Proposition 5.2. For any n G N the law Z„ is the one of a SBM _BtT'^" with drift /i„ conditionally on 
Bj^ .Mil ^ with 

hn{y) = C„exp(B„(y) - S„(a;))p"'^"(r, a;, y), 
where C„ is f/ie normalizing constant such that J hn{y)dy — 1. 
Proof. See [3]. 

Let riQ be fixed and = to < ti < . . . < tng < T. Set yo = 2; to simplify the notations, the law of 
(wfi , . . . , wt„^ , wt) under Z„ is given by 

riQ-l ^ 

hn{v) q'^-^'^{ti+i~t^,T ~ti,yi,y,yi+i)dyi...dynody. (23) 

Once has been sampled along hn{y)dy, we can sample along q^'^" (ii, T, x, wt, yi)(i2/i and each 
wt,,^! along g"'^"(tj+i - ti,T - tj,wt,,WT,J/i+i), using the Auxiliary Algorithm 1. 
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In order to sample along hn{y)dy we make use of the following considerations. We have 

K{,y) =Cn exp (B„(y) - (T, x, y)vi''^" (T, x, y) 

=C„ exp ( - ^i„iy ~^)+ £ Kz)dz) x exp ( + ^l^{y ~ x) ^ ^T)p'''^{T, x, y)v^^^'- (T, x, y) 

=C„e-4^ exp (S(y) - B{x))v^^^- {T, x, y)p'''\T, x, y). 

Let us denote by M an upper bound for the function z i— > ^\{z) (see our assumptions in § 13. ip . Then, 
using the result of Lemma |4. II and performing easy computations, we easily see that for any < (5 < 1 : 

h (it\ „2 1 .1 



with 

Using one may easily check that fg'"''^"{y) < 1 for any y e R. One might then optimize w.r.t. 
S E (0, 1) in order to find /]''"'^" closest to 1. 

Let us set for simplicity, /''■■iT''^" = /1/2 We deduce therefore the following procedure in order to 
sample along hn{y)dy. 



Auxiliary Algorithm 2: Sampling along hn{y)dy 

1. Sample Y ^Af{x,2T). 

2. Evaluate 

3. Draw U ^ U{[0,1]). If U < accept the proposed value Y. Else return to 
Step 1. 



6. Direct exact sampling of a skeleton under Z (Step 2 of the Exact Simulation Algorithm) 



Proposition 12.11 will now play a crucial role. 
Recall the definition ([T]) of 6. Let us denote 



14 



P. Etore and M. Martinez 



1 if 61 > 1 if 6* > 

and Ct 5^ ^ = < 

-f\tM) if ^?<0, \-i\tM)l\T-t,\h\) if d<(\ 

Remember our definitions ([lO]) , (Il2|) , ([13]) and ((17]) and Remark l4?2l It is clear from ([20]) that ^fin ^ 9 (as n 
tends to +00), so that we have, 

vi^''"{t,x,y) > v^{t,x,y) V(i, x, y) e M+ x R x M, 

n— j-oo 

7^'^"(t,z) >l%t,z) V(t,z)eM+xR, 

n— ^-oo 

Ct-r V(t,x)eR+ xM, 



'^r.T.a.b ^ C'f T a 6 ^(^. T, fl, 6) G R+ X R+ X K X 



Let us now examine the sequence of the rejection functions used in the Auxihary Algorithm 2. 

From the same reasons as above, it is clear that (/]''" '^") converges towards 

f, {y) = VT^ exp (^Biy) - - j pO,^T/il - S),x,y) " ^' 

this convergence being dominated. Thus, applying the result of Proposition 12.11 the sequence of laws 
(hn{y)dy) converges to some limit law hg{y)dy. 

In the same manner, for any fixed a, 6 G M, the sequence (/^ ^ " ' " ) of rejection functions used in Auxiliary 
Algorithm 1 converges towards 

.OJ.e/ X v%t,a,y)v'^{T -t,y,b) 
fa,b iy) = TTe < 1> 

this convergence being dominated. 

Consequently, the law q^^f^" (^t, T, a, b, y)dy converges towards a limit law q^{t, T, a, b, y)dy. 

Let again uq be fixed and < ti < . . . < i„Q < T. Passing to the limit in (1^51) we get that the law of 
[ujti , ■ ■ ■ , Wt„j, , i^t) under Z„ converges (with ya — x) towards 

"0 — 1 

he{y) Y{ l^ii^+i - i-uT ~ U,y^,y,yi+i)dyi . . .dy^gdy. (24) 

Consequently, from Proposition 15. 1[ we conclude that the law given by is nothing else than the law 
of (wfi , . . . , iot^^ , lot) under Z. 

Using again Proposition 12 . II and the above considerations we can propose the expected algorithm in order 
to sample skeletons under Z. It will use the two following Limit Auxihary Algorithms. 
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Limit Auxiliary Algorithm 1: Sampling along he{y)dy 



1. Sample Y r~.M{x,2T). 

2. Evaluate 

f^'^Y) < 1. 

3. Draw U ~ W([0, 1]). If U < /^'"(Y) accept the proposed value Y. Else return to Step 1. 



Limit Auxiliary Algorithm 2: Sampling along q^{t,T,a,b,y)dy 



1. Sample a Brownian bridge Y along q°'°{t,T,a,b,y). 

2. Evaluate 

3. Draw U ~ l({[0, 1]). If U < f^f{Y) accept the proposed value Y. Else return to Step 2. 



Performing Step 2 of the Exact Simulation Algorithm. 
Sampling (wti, • . • ,wt„^,WT) under Z (starting from x) 



1. Sample lot along h0{y)dy using the Limit Auxiliary Algorithm 1. 

2. Sample along q" {ti,T,x,(jjT,y)dy using the Limit Auxiliary Algorithm 2. 

3. For i = 2,. . . ,no, sample coti+i along q^{ti+i-ti,T-ti,iJti,iOT,y)dy using the Limit Auxiliary 
Algorithm 2. 
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-6 -4 -2 2 4 6 

Figure 1: Brownian motion with two-valued drift, case 9o = —9i = 1 (T = 1). 

7. Numerical Experiments 

7.1. Exact simulation of a Brownian motion with two- valued (or alternate) drift 

In this paragraph, we choose to exhibit numerical results obtained with the exact limit algorithm for the 
simplest non-trivial cases 

dXt = dWt ± sgniXt)dt, Xo = 0, 

corresponding to either Oq = —9i ~ ±1 in Q (6(y) = ±sgn(y) in ([T])). Indeed, in this symmetric case a 
benchmark is provided by the explicit and computable density of Xt given in [lOJ p. 440-441. 

We draw the renormalized histogram of 10^ samples of Xt and compare it to the explicit density of Xt 
(Figure [1] for the outgoing case 6*0 = 1 and Figure [2] for the incoming case 0o = ^1)- 

In the non-symmetric case we can still use our limit algorithm but the density of Xt becomes less explicit 
(see formula (6.5.12) in |10|). Thus we will use as a benchmark the renormalized histogram of 10^ samples 
oiX^, where {X^) denotes an Euler Scheme with time step A = T.IO"^. We chose = 2, = -1, T = 1 
and X — 0.0. We plot the corresponding renormalized histograms on Figure [3] 
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Figure 2: Brownian motion with two- valued drift, case = —Si ~ (T = 1) 
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Figure 3: Limit algorithm v.s. Euler Scheme for Brownian motion with two- valued drift with 6o = 2 and 8i = —1 
{x = 0.0 and T = 1). 
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Figure 4: Limit algorithm v.s. Euler Scheme for the case where b is given by (|25p {x — 0.0 and T = 1). 

7.2. Exact simulation of an SDE with a discontinuous drift coefficient 

We consider now the SDE ([T]) with 



if a; > 



b{x) 



(25) 



^ - ^ cos 



(fa;) if a; < 0. 



Let < T < oo. We wish to sample along Xt- 
We have 9 = -37r/4 and 



TT 

20' 



We take K — + as an upper bound for 0. This allows to use the limit Algorithm. 

Figure [3] shows a comparison between a renormalized histogram of 10^ samples of Xt obtained with the 
exact limit algorithm, and a renormahzed histogram of 10^ samples of X^ , where (X'^) denotes an Euler 
Scheme with time step A. We chose x = 0.0, T 1 and time-steps A = T.IO"^ and A = T.IO"'". 
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